Test of cold denaturation mechanism for proteins as a function of water's structure 
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In a recent paper [PRL 91, 138103 (2003)] a new mechanism to explain the cold denaturation of 
proteins, based on the loss of local low-density water structure, has been proposed. In the present 
paper this mechanism is tested by means of full atom numerical simulations. In good agreement 
with this proposal, cold denaturation resulting in the unfolded state was found at the High Density 
Liquid (HDL) state of water, at which the amount of open tetragonal hydrogen bonds decreases at 
cooling. 



Hydrophobicity is an important driving force in the 
process of polymer folding and protein stabilization 
E II II 11 II II 0, II E3 Basically, hydrophobicity pro- 
vides a mechanism to avoid the energy cost of disruption 
of a water hydrogen bond by the presence of an apolar 
molecule . Actually, the folding of a polymer is partially 
due to the assembly of the apolar monomers in order to 
minimize the disruption of hydrogen bonds among water 
molecules ^HtHl- Also, the stability of the folded state 
is enhanced by the cage structures formed of highly orga- 
nized water molecules around the biopolymer where the 
disruption of hydrogen bonds is minimized . Increas- 
ing the temperature of a polymer in a water bath leads 
to structures different from the native folded state. Hy- 
drophobicity is reduced by thermal fluctuations and the 
proteins stability is lost. However it is also possible to 
reduce polymers hydrophobicity and to induce unfolded 
states by lowering the temperature, a phenomenon called 
cold denaturation 0, 0, ^| . The microscopic mecha- 
nism leading to cancellation of the hydrophobic effects 
and cold denaturation is not clear, due in part to the 
complexity of the protein-solvent interactions. Cold de- 
naturation has two important properties, (i) For pure 
water solvent, cold denaturation is only found at very 
high pressures (P > 2kbar) So, in order to get high 
pressure cold destabilization of the folded protein, hy- 
drophobicity should decrease when cooling at high pres- 
sures, (ii) The transition to the unfolded state is driven 
by a smaller volume, indicating a positive slope in the 
P-T transition line. On the contrary, the usual thermal 
denaturation behaves with negative slope. Hydrophobic- 
ity weakening at high pressures should also be accompa- 
nied by a positive slope. In recent years there has been a 
huge amount of work in order the study the dependence 
of hydrophobicity with pressure. Contact and solvent 
separated configurations have been studied from an en- 
thalpy/entropy point of view Solvent separated con- 
figurations are increasingly stabilized at higher pressures 
by enthalpic contributions. Full atom simulations of two 
Lennard- Jones particles in model water have found that 
the aggregation of the two particles becomes unstable 
for pressures in the kbar range 0, and recent replica 



exchange molecular dynamics simulations have studied 
the weakening of the hydrophobic interaction with pres- 
sure for a C-terminal fragment of protein G These 
findings are in agreement with a cancellation of the hy- 
drophobic force at high pressures. Actually some recent 
molecular dynamics simulations have already found the 
swelling of a protein when lowering the temperature at 
high pressures with a positive transition slope as corre- 
sponding to the cold dentauration of real proteins [l^ . 
However the question still remains: Why is hydrophobic- 
ity canceled at very high pressures in such a way as to 
promote cold denaturation? 

During the last years, several works have attempted to 
answer this question by using different models: Thermo- 
dynamical models [2£j, effective attraction models |2lj . 
temperature dependent attraction models 22], effective 
water-protein interactions models including solvent ex- 
clusion 0, El EH ? lattice models |2^ based on a bimodal 
description of the energy of water in the shell around the 
hydrophobic molecule |2£( and models mimicking the in- 
teraction between water and non polar monomers |2f| . 
Very recently, a lattice model that captures the effect 
of volume correlations on water at high and low pres- 
sures has been used to propose a possible mechanism 
for cold denaturation and hydrophobicity weakening at 
high pressures [2jj. The proposal is based on the known 
existence of two possible water metastable states. A 
low density /pressure liquid with a local structure very 
similar to the low-density amorphous (LDA) solid wa- 
ter, consisting of an open tetrahedral network, and a 
high density /pressure water, very similar to high-density 
amorphous (HDA) solid water, where core-core repulsion 
dominates, similar to simple, non anomalous liquids |2^ . 
A way to change from the low density liquid (LDL) to 
the high density liquid (HDL) is by increasing the pres- 
sure at constant temperature. Actually, there is a critical 
pressure P c separating both behaviors. For P < P c a de- 
crease in temperature implies an increase in the number 
of energetically favorable open tetrahedral interactions, 
and an increase on volume. For a given temperature 
(called temperature of maximum density (TMD)) this 
increase in volume is more important than the regular 
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compresion of the liquid, and the density reaches a max- 
imum. Then, the density decreases more and more until 
reaching the negative sloped Ice Ih crystalline line ,or 
the negative sloped Kauzmann boundary (vanishing dif- 
fusivity) in numerical simulations. On the contrary for 
P > P c the number of interactions forming open tetrahe- 
dral structures decreases when cooling down the system, 
producing the more dense interpenetrating tetrahedral 
network and the more energetically dominant core-core 
interactions. Now there is no TMD, and the density al- 
ways decreases (similar to a non-anomalous liquid) until 
reaching the positive sloped, Ice III or Ice V crystalline 
line (or the positive sloped Kauzmann boundary in nu- 
merical simulations). Basically, for a constant tempera- 
ture, the higher the pressure, the lower the number of low 
density open tetrahedral structures. The hypothesis pro- 
posed in ref. |2jj is that hydrophobicity (i.e. aggregation 
to avoid the disruption of hydrogen bonds and formation 
of open cage structures) is drastically reduced for the 
HDL state. Basically, at the HDL state, the existence 
of solvent-separated apolar molecules is not energetically 
unfavorable as core-repulsion interaction is now the pre- 
vailing force instead of the open tetrahedral hydrogen 
bonded network. Also, the stability due to the formation 
of low density hydrogen bonded cage structures is not 
energetically favorable anymore. The high density wa- 
ter molecules are now capable of penetrating inside the 
protein leading to unfolding With this hypothesis is 
easy to understand why cold dentauration in pure water 
is found at P > 2kbar, since P c ~ 2kbar in pure water. 
Also it is possible to understand why the cold denatu- 
ration slope at the P-T phase diagram is positive, since 
the non-anomalous, high pressure, ice crystalline lines 
(such as Ice III or Ice V) have also a positive slope (at 
P > P c the density of water always decreases when cool- 
ing down). In order to check this idea, a simple lattice 
model was proposed in Ref. [2jj and the results were suc- 
cessfully compared with the ones corresponding to real 
proteins such as Ribonuclease A (2^. However, a direct 
comprobation of this hypothesis by means of "real" full 
atom simulations was still lacking. The aim of this pa- 
per is to present results from full atom simulations which 
clearly correlate the weakening of hydrophobicity to the 
presence of the high density state of water. 

In this paper the SPC/E model of water is simu- 
lated at different densities and temperatures. For posi- 
tive pressures the line of temperatures of maximum den- 
sity (TMD) has a negative slope in the P-T plane and 
disappears for pressures P > P c ~ 300MPa |3l| . From 
the point of view of density-temperature simulations it 
implies that there is a critical density p c ~ l.lg/cm 3 . 
If we do constant density simulations at p < p c we 
are at the LDL regime, and the number of open hy- 
drogen bonded tetrahedral networks grows as the tem- 
perature lowers. On the contrary, if we are at p > p c 
we are at the HDL regime and the number of open hy- 
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FIG. 1: Density-temperature phase diagram of SPC/E water. 
Black points correspond to the TMD line and white points 
correspond the the Kauzmann boundary 31]. Large crosses 
indicate the p vs. T conditions considered. Note how both 
regimens (LDL and HDL) are taken into account. 



drogen bonded tetrahedral networks decreases as low- 
ering the temperature. For this reason , in order to 
study both regimes, we perform numerical simulations 
at p = lg/cm 3 ,l.lg/cm 3 ,1.2g/cm 3 and 1.3g/cm3 (see 
Fig.l). At all densities, we perform simulations at T = 
175K, 225K, 275K and 325K in order to sample from 
highly thermally activated water molecules to conditions 
close to the Kauzmann boundary. First, in order to con- 
firm the correctness of the density-temperatures phase 
points chosen for the simulations, the oxygen-oxygen pair 
correlation function for a sample with a maximum num- 
ber of 460 water molecules has been measured. A re- 
action field technique 32] with a cut-off of 0.79nm ac- 
counts for the long-range Coulombic interactions. Con- 
figurations were sampled by means of a non-constant 
step Monte Carlo algorithm. Results are shown in Fig. 2. 
They match previously reported calculations and re- 
cent experimental results 33]. At p < p c a maximum 
of the oxygen-oxygen correlation function at r ~ 4. 5 A is 
found. This maximum grows as we cool down the sys- 
tem and is typical of LDA water and open tetrahedral 
network structures (like Icelh). It has been also mea- 
sured by numerical simulations in LDL water 28]. On 
the contrary, at p > p c , the peak at r ~ 4. 5 A decreases 
when cooling down the system an a new peak at r ~ 3. 5 A 
rises up. This new peak is typical of the HDA water |3^ 
and interpenetrating dense tetrahedral networks (like Ice 
VI and Ice VII). It has also been measured by computer 
simulations in model water [2^. With this checking we 
are now confident about the range of pressures to be sim- 
ulated, and we know when we are dealing with LDL and 
when we are dealing with HDL. 
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FIG. 2: Oxygen- Oxygen pair correlation functions vs. dis- 
tance (in Angstroms) for the p vs. T conditions simulated in 
the paper. Thicker lines correspond to lower temperatures. 
Vertical lines indicate the values of r were main secondary 
peaks appear: r ~ 3.5A for the HDL and r ~ 4.5A for the 
LDL. 



To test the hypothesis proposed in Ref. |2jJ we should 
now embed an apolar polymer chain in our SPC/E wa- 
ter bath. If the hypothesis is correct, at low temper- 
atures, folded states are expected for p < p c and un- 
folded states are expected for p > p c . In order to at- 
tribute any configuration change to the influence of the 
solvent, we simulate a modified model of the polyethy- 
lene molecule 34], free of bond angles and dihedral po- 
tential barriers. The polymer-chain consist of N = 15 
hydrophobic CH2 beads, represented by Lennard- Jones 
interaction sites with e = 0.112kcal/mol and a = 1.33 A 
and linked to quasi-rigid harmonic bonds of 1.53 A length 
and interaction coupling k r = 350kcal/mol. The water- 
polymer cross parameters were obtained using the con- 
ventional Lorentz-Berthelot mixing rules. The structure 
of the hydrophobic polyethylene chain is characterized by 
its radius of gyration Rq = 1/N^2((fl — rw) 2 ) (being 
r^c = l/N^fi). For the extended polymer the radius 
of gyration is Rq ~ 4A and a typical folded configura- 
tion with spherical-like symmetry has a gyration radius 
close to Rq = 2.75A. So, we consider a configuration 
to be folded if the gyration radius is Rq < Ro = 2. 75 A 
(see Fig. 3). For each pressure- temperature point of the 
phase diagram we measured up to 10 8 thermally equi- 
librated polymer configurations in order to build a his- 
togram with the most probable values for the gyration 
radius. Using these data we are able to calculate an av- 
eraged value and the corresponding dispersion. Results 
are shown in Fig. 4 (dispersions are represented as error 
bars). 

Clearly there are two different behaviors. For p < p c 




FIG. 3: (a) (Color online) Typical folded configuration of 
the polyethylene chain Rq ~ 2A and (b) typical unfolded 
configuration Rg ~ 3.5 A. Atom's sizes are not at scale. 



and low temperatures (T < 275K) the protein adopts 
a folded state with a gyration radius Rq < Ro (see 
Fig. 4). There are no significant changes by cooling down 
the system. The only way to obtain unfolded states is 
by increasing the temperature up to T = 325K. In 
this range of densities of the system, corresponding to 
the LDL state, no cold denaturation is found. On the 
contrary when we increase the pressure up to p > p c 
(HDL state) there is a dramatic change on the protein 
status by decreasing the temperature (see Fig 4). At 
low temperatures, the gyration radius grows to a value 
corresponding to the unfolded state. We may defini- 
tively say that, in this range of densities of water cor- 
responding to HDL, cold denaturation is found at low 
temperatures. In order to get a more direct relation be- 
tween the value of the gyration radius and the water's 
tetragonal net status, it is possible to define a parame- 
ter (S = g - (r = 4.5A) - g _ (r = 3.5A)/(g _ (r = 
4. 5 A) + g - (r — 3.5A))) based on the two secondary 
maxima of the oxygen-oxygen pair correlation function. 
If the water's status corresponds to a LDL, the value of 
the correlation function at r = 4. 5 A is going to be much 
important than the value of the correlation pair distri- 
bution function at r = 3. 5 A (see Fig. 2), and S is going 
to be positive. If the status corresponds to a HDL, the 
value at r = 3. 5 A is going to be larger than the value at 
r = 4.5A (see Fig. 2) and S turns negative. Fig. 5 is a 
plot of the different gyration radius values found Rq vs. 
S (see Fig. 5). The set gets divided in four quadrants: (i) 
Rg > Ro and S > (unfolded state and LDL water), (ii) 
Rg > Ro and S < (unfolded state and HDL water), 
(in) Rg < Ro and S > (folded state and LDL water), 
(iv) (1) Rg < Ro and S < (folded state and HDL 
water). (In order to avoid effects on the water structure 
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FIG. 4: Mean value and dispersion of the gyration radius (in 
Angstroms) for the polyethylene chain at the different p vs. 
T conditions considered (dotted line is a guide to the eye). 
Continuous line indicates the value Ro = 2. 75 A. At p = 
lg/cm 3 , l.lg/cm 3 (LDL state) no cold denaturation is found 
by cooling and the gyration radius remains almost constant. 
At p = 1.2g/cra 3 , 1.3g/cm 3 (HDL state) cold denaturation is 
found and the gyration radius grows when cooling down. 



due to thermal fluctuations we focus our attention on the 
low temperature regime T < 275K). Note how almost 
all data corresponds to cases (ii) and (in ) , in agreement 
with the hypothesis proposed in Ref. |27j . 

To conclude, full atom simulations of the folding be- 
havior of a model polyethylene embedded in SPC/E wa- 
ter show that there is a correlation between the cold de- 
naturation and the density of water. For LDL the exis- 
tence of a local structure consisting of an open tetrahe- 
dral hydrogen bonded network prevents cold denatura- 
tion. The volume of this network increases when cooling 
down the system and so does hydrophobicity. On the 
contrary, for the high pressure HDL, the local structure 
is a packed tetrahedral hydrogen bonded network dom- 
inated by core-core repulsions. The volume of this net- 
work decreases when cooling down the system and so does 
hydrophobicity, explaining the existence of high pressure 
cold denaturation of proteins. 

Comments from J.A.Gonzalo and H.E.Stanley are 
greatfully acknowledged. 
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